Stability of excited states of a Bose-Einstein condensate in an anharmonic trap 
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I. INTRODUCTION 



The mean field theory of a Bose-Einstein condensate 
(BEC) is based on the Gross-Pitaevskii equation (GPE) 

1 



i% = -A* + V(x)$ - ct|*| 2 *, 



(1) 



for the macroscopic wave function 'J = x). This 
model describes accurately the behavior of an ultracold 
condensed atomic cloud trapped by an external poten- 
tial V(x). In the dimensionless variables of Eq. (JXJ) the 
Planck constant is equal to one, h = 1, and the atomic 
mass is m = 1/2. The parameter a stands for the sign op- 
posite to that of the scattering length a s : a — —sign (a s ). 
An important class of solutions of GPE are stationary 
nonlinear modes defined as 



*(t,x) = e- iw V(x), 
with the boundary conditions 

V>( x ) ~> as l x l — * °°- 
In this case i/j(x) satisfies the equation 

Aip + (ui- V(x))ip + ctV> 3 = 0. 



(2) 
(3) 
(4) 



state solutions of Eqs. ©-(H]) (i.e. its positive solutions 
which minimize the energy functional for Eq. ([1])) are 
of primary importance for BEC applications [ij . Apart 
from them some non-ground nonlinear modes have also 
been studied (see e.g. [1 [1 0, 0, & EH)- 

However all of 

the practical applications of high order modes are linked 
to their experimental feasibility, what requires stability. 
The stability of high order modes has already been stud- 
ied in the case of a harmonic potential F(x) = |x| 2 . The 
one-dimensional case was studied in Refs. [H EH El EH 
while multidimensional solutions were considered in Refs. 

[H [M m El El Hi EE H HI • 

However, to the best of our knowledge, the relation 
between the stability properties of nonlinear modes and 
the specific forms of the confining potential V(x) has not 
been discussed, previously. This is a problem of signifi- 
cant practical importance since high order modes can be 
used for generation of nonlinear coherent structures, such 
as for example solitonic trains in quasi-one-dimensional 
limit (see e.g. [Toj). 

In this paper we study how the shape of the poten- 
tial V(x) governs the stability properties of nonlinear 
modes. In our study and to simplify the analysis we 
will consider a quasi-one dimensional geometry modeled 
by a one-dimensional GPE [14| 



The nonlinear eigenvalue uj is called the chemical po- 
tential in the context of BEC applications. The ground 



In that situation, Eq. (01 becomes 

ip xx + (w - V{x))i> + enjj 3 = 0. 



(5) 



(6) 
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We will show that the harmonic potential V(x) = x 2 
corresponds to a very peculiar situation related to the 
fact that in the linear limit this potential has an equidis- 
tant spectrum. We will discuss how switching from the 
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harmonic potential to an anharmonic one makes higher 
nonlinear modes "more stable" and even a "weak" an- 
harmonicity (say, V(x) = x 2 + nx 4 , \n\ <C 1) is enough to 
change drastically the stability properties of high-order 
nonlinear modes. 

It is relevant to point out that the nonlinearity intro- 
duces not only mathematical difficulties for the study 
of the eigenmode problem, but significantly diversifies 
the list of physically relevant limiting cases. Indeed, the 
linear case has two characteristic scales: the de Broglie 
wavelength and the scale of the potential. In our no- 
tations, the first one is A ~ max x \tpx\/ max x \ip\; the 
scale of the potential, roughly speaking coincides with 
the classically allowed domain (in what follows denoted 
as Xd)- The nonlinearity introduces a third relevant scale 
I ~ l/max^l^l- For repulsive nonlinearities this scale 
corresponds to the healing length [l), while in the case 
of attractive nonlinearities it measures the width of a 
matter-wave soliton. Therefore, the diversity of limiting 
cases of the nonlinear eigenvalue problem is characterized 
by the interplay between the parameters A, Xd, and I. 

The paper is organized as follows. In the introductory 
section |TT] we describe the typical structure of the families 
of nonlinear modes and formulate the stability problem. 
The results for the harmonic potential are summarized 
in Sec. IIII1 Next, we turn to the consideration of the 
anharmonic potentials V(x) = x 2 + nx 4 (Sec. IIV[) and 
V{x) = x 4 (Sec. [V} . The last section PVTl contains some 
additional discussions and a summary of our results. 

II. DEFINITIONS AND PREVIOUS RESULTS 

A. Branches of nonlinear modes 

In the small amplitude limit tp — > 0, the cubic term 
in Eq. © can be neglected and the solutions can be 
approximated by the eigenfunctions ip n (x), n = 0, 1, . . . 
of linear eigenvalue problem 

i> xx + (w - V(x))i> = 0. (7) 

It is assumed that the potential V(x) is nonsingular, 
bounded from below, and V(x) — > oo as \x\ — > oo, so 
that the spectrum of (J7J is discrete. Throughout this 
paper we will deal with even potentials, V(x) — V(—x). 
The real eigenfunctions of (0, we denote them ip n (x) 
with n — 0, 1, . . ., constitute an orthonormal set 

/oo 
■4>n(x)lpm(x) dx = S m , n (8) 
-oo 

(here S mjn is the Kronecker delta). 

It will be convenient to describe the families of nonlin- 
ear modes in terms of bifurcation diagrams in the plane 
(u>,N) where 

/oo 
^ 2 (x)dx, (9) 
-oo 



corresponds to the number of particles. Then, the re- 
spective eigenvalues of (J7J, u = u> n , n = 0,1,..., are the 
points of bifurcation where families of nonlinear modes 
of Eq. ©, to be denoted as To, Tx, T2, . . ., branch off 
from the zero solution ^(x) = 0. The branching off takes 
place for both cases, a — 1 (attractive nonlinearity) and 
(j = — 1 (repulsive nonlinearity). Therefore we will label 
the branches of nonlinear modes in the corresponding 
cases by and Tn\ At the same time, in the state- 
ments which are applicable to both cases a — ±1, we 
omit the superscript writing simply r„. Two examples 
of diagrams for V(x) = x 2 are shown in Fig. [TJ Following 
Q we call the modes described above modes with linear 
counterpart, since they can be viewed as modes of a linear 
oscillator "deformed" by the action of the nonlinearity. 

A simple analysis shows that in the vicinity of a bifur- 
cation point, say u — uj n , the small-amplitude solution 
of Eq. for the branch T n can be described by the 
asymptotic expansions 

ip n (x) = eip n (x) + o(e), (10a) 
w n = Q n - e 2 aVt n + o(e 2 ), (10b) 

where e 1 is a small parameter and the coefficient Q n 
is given by 

/oo 
$*{x)dx. (11) 
-oo 

From the physical point of view, Q„ describes the two- 
body interactions, and thus defines the characteristic 
scale I. Hence the limit described by Eq. (fT0|) in physi- 
cal terms can be defined as I <C A. We observe that this 
limit can be achieved not only due to small numbers of 
condensed particles, expressed by the condition e <C 1 
[due to the normalization ((SJ)], but also for large enough 
n, since fl n — > as n g rows. For instance, for the case 
V(x) = x 2 we have [31( 

Q n = 0(n- 1 / 12 ). (12) 

B. Physical units 

Throughout this paper the nonlinear modes will be 
characterized mainly in terms of the number of particles 
N and the mode frequencies lu. Since we will be inter- 
ested in applications of our results to the BEC mean 
field theory while our analysis will be carried out in di- 
mensionless units, before going into details we outline 
the link of our variables with experimentally feasible pa- 
rameters. To this end we notice that the dimensionless 
form of Eq. © corresponds to the situation where the 
distance and time are measured respectively in units of 
ao and 2/u)q, where ao is the longitudinal length of the 
trap and loq = fi/(mag), with m being the atomic mass, 
is the " effective" longitudinal trap frequency (in the case 
of a parabolic potential it is the real frequency of trap, 
while oq is the linear oscillator length). In the chosen 
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scaling the energy is measured in the units huJo/2. Then 
it is a straightforward algebra to ensure that the link be- 
tween the real, i.e physical, number of particles Af and 
the norm of the solution N (also referred to as num- 
ber of particles) introduced by ([9]) is given by the for- 
mula N = (a±/A^/2Tr\a s \)N, where a± is the transverse 
linear oscillator lengths and a s is the s-wave scattering 
length. In this last formula, as well as in the reduc- 
tion of the 3D Gross-Pitaevskii equation to the ID model 
we have assumed that the trap is cigar-shaped, i.e. 
a± <C do, and that in the transverse direction the trap 
is harmonic. Now it is not difficult to estimate, that, for 
example for a trap with characteristic values ao ~ 1 mm, 
aj_ ~ 10 /im, and a s ~ 5 nm, a unit of the norm ./V corre- 
sponds to TV ~ 10 2 atoms, or to the mean atomic density 
- 0.4- 10 9 cm- 3 . 



C. Quasiclassical quantization 

The limit of large N and small local densities \ip\ 2 -C 1, 
denoted nonlinear WKB approximation 0], allows for 
an analytical construction of explicit solutions. Let 
us consider V(x) — x 2d , where d is a positive inte- 
ger and assume that N 3> 1, or more precisely, that 



8 = N^ 1 ^ 1 ^ <C 1. Bound states corresponding to suffi- 
ciently large u, specifically E — u)N~ 2 > 1, correspond 
to an atomic cloud distributed over a large spatial do- 
main roughly determined by the classical turning points, 
±Xd = ±uj 1 / 2d . Since x d grows with uj, one can reach 
the levels corresponding to sufficiently low densities of 
particles, i.e. the quasi-linear limit. 

Returning to the definition of modes with linear coun- 
terpart, the arguments presented in this section allow one 
to conjecture, that only modes with a linear counterpart 
can exist in the limit uj — ► oo at N fixed. 

To describe that situation, we focus on the repulsive 
case a < 0, introduce a new independent variable £ = 
g 1 /( 1 + d ) x anc j a renormalized macroscopic wave function 
4>{Q = <S (d - 1)/2(ci+1) VK£), and rewrite Eq. © as follows 

5 2 cj> C( + (E- ( 2d )0 - 6<P 3 = 0. (13) 

Since now 6 *C 1 Eq. (|13[) is a convenient represen- 
tation of the stationary eigenvalue problem for the ap- 
plication of the nonlinear WKB approximation. Skip- 
ping details, which can be found in [9j, here we present 
the equation implicitly defining the diagram in the plane 

(E,sy. 



E ^ K A d Y\ n+ V + 7B hx {— )-° dE 2d+&n 

Here n stands for the energy level and we introduced the constants 



1 - i 



B rl 



(14) 



A d 



1 n-l 

yl - y 2d dy 

-l 
tZ-l 



^ = 35^^008 ^111, 



fc=l 



B *=2d 



1 - cos [if) 
1 + cos (if) 



dy 



-i y/i - v 2d 
kit 

11 ' ~d~ arc ^ an 



for d= 1,2,.. 



ford = 2, 3,... 



and C\ = 0. T(-) is the standard notation for the gamma 
function 1 301. 



the nonlinear WKB equation (fT4|) acquires the form 



£ 3 / 4 = 



Inn 



2V2nK 



tt \ n + — 



O 



(15) 



Solutions of the transcendental equation ([14]) with re- 
spect to the energy E at fixed S and n give the eigenval- 
ues (energy levels). According to the previous discussion, 
when E — > oo one recovers the WKB formula for the en- 
ergy levels of the potential V(x) — x 2d . While for the 



case of the harmonic oscillator the form of ([14)) can be 
found in Ref. Q , for the situation of our particular inter- 
est below, d — 2, the limit E — > oo (and thus n — > oo) of 



where K(-) is the complete elliptic integral of the first 
kind [U. 

Considering the last two terms in the expansion of the 
energy levels (I15[) as a perturbation for n large enough, 
and neglecting them in the leading order, in remaining 
part of the expression for E n one readily recognizes the 
familiar Bohr-Sommerfeld quantization condition. Thus 
formula (|15p can be viewed as the nonlinear general- 
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ization of the standard quasi-classical quantization well 
known in the quantum mechanics. 

D. The stability problem (general) 

Consider now the stability problem for the modes cor- 
responding to a fixed branch r„. Let ipn(x) be a solution 
of Eq. ([6|) corresponding to w n . Following the standard 
procedure we represent ^(x, t) = (ip n {x) + £_(x,t))e~ lUnt , 
linearize the dynamical equation with respect to £(x,£), 
and arrive at the equation 

iff = -£«» - (w„ - V(x))£ - a^ n (2£ + £*) , 

where the asterisk stands for the complex conjuga- 
tion. Decomposing into real and imaginary parts, 
£(x, t) = x( x > t) + i<p{x, t) we obtain 

Xt = ~L~ <p, (ft = L+ x, (16) 

where 

L+ = L„ + 2,ail) 2 n (x) , L- = L n + aip* (x) , 

and 

d 2 

L n = ^ + w " - v ( x )- 

It follows from (|16p that stability of the nonlinear mode 
is determined by the spectrum of the following eigenvalue 
problem 

L-L+ ( = A C- (17) 

Since the operator L~ is degenerate (the kernel of this 
operator contains at least the function ^j n (x)), the eigen- 
value problem (|17p has a zero eigenvalue. Then if the re- 
mainder of the spectrum of L~L^ is real and nonnegative 
then the nonlinear mode ^(x, t) = e~ luJnt ip n (x) is said to 
pass the linear stability test. The presence of negative 
or complex eigenvalues in the spectrum of L~L^ implies 
the linear instability of this mode. 

E. The stability problem (small amplitude modes) 



The operator C n is self-adjoint and its spectrum consists 
of eigenvalues \ n ,k = w n — k = 0, 1,2, . . .. In the 
limit e = one has the operator L~L^ — L 2 n and its 
spectrum becomes A n ^ — (cD„ — ^fc) 2 , k = 0,1,2,.... 
Since all A n ^ are real and nonnegative, then if there 
are no multiple eigenvalues in this spectrum, the small 
amplitude nonlinear modes are linearly stable for both, 
repulsive and attractive nonlinearities. However, if the 
spectrum of C„ includes multiple eigenvalues the stability 
analysis implies the study of splitting of these ei genv alues 
when passing from £ = 0to0<e<§;l (see e.g. [H, Hj)- 

F. Krein signature 

Let n be fixed and a pair (A, £(x)) be a solution of 
eigenvalue problem (|17|) where A > is a semi-simple 
eigenvalue of L~L^ and the corresponding eigenfunction 
C(x) is real. It is useful to assign to any such a pair 
(A, ((x)) the value 

Assign (L+C(x),C(a;)) 

called the Krein signature (25|. As the solution ip n (x) of 
Eq. {BJ varies along the branch r„ together with u> n , the 
eigenvalues of the operator L^L~ also vary, but the Krein 
signature of any pair (A, is conserved while there 

is no collision between eigenvalues. When a collision 
between a pair of real positive eigenvalues takes place, 
they can become complex only in the case when their 
Krein signatures are opposite; otherwise these eigenval- 
ues pass through each other both remaining real. So, as 
io n varies along the branch T n the interactions of eigen- 
values with opposite Krein signatures may affect the sta- 
bility of modes in this branch. 

The inverse statement is also valid but in a generic 
situation only [25]: if the Krein signatures of colliding 
eigenvalues are opposite then generically after collision 
they become complex. However additional symmetries 
of the solution can destroy this picture: it will be shown 
that in some cases eigenvalues with opposite Krein signa- 
tures can also pass through each other without causing 
instabilities. 



If lu lies close to a bifurcation point lu ti then the spec- 
trum of the operator L~L~^ can be analyzed by means of 
the asymptotic expansions p0|. Specifically, 

L+ = C n + ae 2 (3^ 2 n (x) - Q n ) + o(e 2 ), (18) 

L~ - C n + o-e 2 {4> 2 n (x) - Q n ) + o(e 2 ), (19) 

L~L+=C* +s 2 aM n + o(e 2 ), (20) 



where 



Cn = dx? +^~V(x), 



(21) 



M„ 



3C n ^ 2 n {x) + (^ 2 n (x)-2n n )C n . (22) 



III. RESULTS FOR THE HARMONIC 
POTENTIAL V(x) = x 2 . 

A. General comments 

In the case of the harmonic oscillator, where V(x) = 
x 2 , the branches r„, n — 0, 1, 2, . . . are depicted in Fig. [1] 
(a), (b). It follows from Fig. Q]that the branches T n arc 
represented by monotonic (at least for moderate values of 
N lu and n) functions N(ui). Previous numerical results 
\V\\ allow to conjecture that there are no solutions with- 
out linear counterpart for this potential. It is known that 
the solutions corresponding to the branches Tq and Ti in 
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FIG. 1: (a,b) The lowest branches of the nonlinear modes of 
Eq. (O for the potential V(x) = x 2 for repulsive (a — —1) (a) 
and attractive (a — 1) (b) nonlinearities respectively. (c,d) 
The same branches but for V{x) — x 4 . All the modes bifur- 
cate from the linear harmonic oscillator modes corresponding 
to the limit TV — > 0. The fragments of curves corresponding 
to stable solutions are shown in bold. 



both, attractive and repulsive cases, are stable (see, e.g. 
for instance [hJ HH and [HJ for more detailed analysis 
of perturbed solutions from Numerical calculations 

show that the modes from (the repulsive case) are 

unstable. In the attractive case the family corre- 
sponds to unstable modes for ix> close to the bifurcation 
point ll>2 — 5. More specifically, in Ref. [ll[ the instabil- 
ity has been observed for uj* < ui < 5 where uj* w 3.83. 
It has also been found that for uo < lu* the mode is stable. 



case e = 0. The spectrum of the operator C n is equidis- 
tant and consists of the eigenvalues A n ,fe = 2(n — k) 
and the corresponding eigenfunctions are given by ipk{x), 
k = 0, 1, . . .. All the eigenvalues are simple and there is 
one zero eigenvalue. The eigenvalues of the operator £^ 



are A r . 



4(fc 



and they correspond to 



the same eigenfunctions tjjkjx), k — 0, 1, . . .. This me ans 
0-th 1-st 2-nd 3-rd 4-th 5-th . . . 
Eigenf-n %j)o{x) i>i(x) 4>2(x) ^(x) ip4(x) i>s(x) ... 



2 
Cl 



2 0-2-4 -6 
|T| [T] 16 36 



64 



£2 4 2 -2 -4 -6 

|~L6~| |T| |T| [w\ 36 



2 
£2 



TABLE I: Eigenvalues of ci,2 and ci,2- 
to emphasize the doubles eigenvalues 



The boxes are used 



that the spectrum of Ct^ includes n double eigenvalues 
A Ui k = 4(rt — k) 2 , k = 0, 1, ... (n — 1), one simple zero 
eigenvalue and infinitely many simple positive eigenval- 
ues. The mechanism of emerging of double eigenvalues 
becomes transparent from Table |TJ Each of the double 
eigenvalues A n ^ has an invariant subspace spanned by 
two functions, ipk( x ) an d ip2n-k(%)- 

Following Ref. [27} we consider the 2x2 matrices 



M 



n,k 



(M n ip k ,'ip k ) (M„^ fe ,-0 2 „-fc) 

(M n 'ip2n-k,'ipk) (Mn^n-^^B-fe) 



B. Small-amplitude modes 

Let us now analyze the stability of the branch T n , for 
both attractive and repulsive cases when uj n lies close to 
the bifurcation point Cj n . In the bifurcation point (linear 
limit) the solutions of the eigenvalue problem Q are the 
pairs (u> n , ij) n (x)), n = 0, 1, . . . where 

u n = 2n + l, (23a) 

Mx) = 1 H n {x)e-i x2 (23b) 
y / 2 n nl y /ir 

and H n {x) is n-th Hermite polynomial (see e.g. [3(|). 

The stability of small amplitude solutions can be stud- 
ied using the formulas l|18 p -(|22 [) . Let us start with the 



If the eigenvalues of M nM are fi^ k and fi^ k , fi^ k ^ fj^ k , 
then for e <C 1 the double eigenvalue A„jt of C? n splits 
into two simple eigenvalues of L~L+: A~\ = A„^ + 
e 2 cr/i~' fe + o(e 2 ) where j — 1,2. Therefore, if the eigen- 
values of any of the matrices M ni k, k = 0, . . . , n — 1 are 
complex, the instability of the small-amplitude solution 
ty(t,x) — e~ luJnt ij) n (x) takes place. It is important that 
since both, repulsive (a = —1) and attractive (a = 1) 
cases are described by the eigenvalues of the same ma- 
trices M ra ,fc, the complex eigenvalues of M„,fc for some k 
means the instability of small-amplitude modes in both, 
repulsive and attractive cases. 

Simple, but tedious algebra gives the expressions for 
the elements of the matrices M n ,k' 
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(n-k) 



Tr2( n+k )nlk\ 



(M n ipk,ipk) 

{M n i) k ,^2n-k) = -(M n 1p 2 n-k,1pk) = 
{M n ^2n-k^2n-k) 



Hl{x)H 2 k {x)e 



dx 



4(n - fc) 
7r2 2n (n!) 2 



H*(x)e 



2x' 



dx, 



(n — k) 



4(n-k) 



7r 2(3«-fc) n !(2n~ fc)! 



Tr2 2n n\y/kl{2n- fc)! 



H l{x)Hl n _ k {x)e 



Hl(x)H 2 n-k{x)H k (x)e 2x dx, 



dx + 



4(n - fc) 
7r2 2 "(n!) 2 



H*(x)e 



-2x s 



dx. 



n = 1 n = 2 7i = 3 n ; 



fc = 



0.199 




fc = 1 



fc = 2 



fc = 3 



fc = 4 



fc = 5 



C 



C 



0.125 




c 



c 



c 



0.089 




c 



c 



0.068 




0.055 




n = 6 



0.133 -0.005 -0.135 
0.404 0.618 0.816 



0.058 
0.473 



C 



C 



0.045 




TABLE II: Double eigenvalues of the matrices M n ^. Each 
cell of the table for n > fc contains either letter "C" , meaning 
that the eigenvalues are complex or two real eigenvalues. 




Using Maple we calculated the eigenvalues of the matrix 
M n The results are collected in Table [TTI were we ob- 
serve the following facts: 

(i) In columns 2 to 6 there is at least one letter "C" 
which means instability of small-amplitude modes be- 
longing to the respective branch r„, for both, attrac- 
tive and repulsive nonlinearities. We conjecture that the 
instability of small-amplitude modes takes place for all 
branches and with n > 2. 

(ii) In the case n = 1 the mode is stable. That con- 
firms the results of [HI (see Fig. 1 there). It is interesting 
that in the limit e = the algebraic and geometric mul- 
tiplicities of the eigenvalue A = 4 are both equal to 2. 
At the same time two real eigenvalues emerging from the 
double eigenvalue A = 4 for e -C 1 have opposite Krein 
signatures which do not correspond to the generic case 
(see 1H). 

(iii) When n = k + 1 the matrix M n ,k has a zero eigen- 
value. This reflects the fact that 



CO*) 



dx 



(24) 



is an eigenfunction of the operator L n corresponding 
to A = 4 for any mode ip n {x) belonging to any branch 



FIG. 2: Branches ri a ' r ' and plots of real and imaginary parts 
of eigenvalues A of the operator L~L^ vs u for V(x) — x 2 and 
n — 2. (a) N vs lo for the branch F^"' (the part of the branch 
corresponding to stable solutions is shown in bold), (b) Real 
and (c) imaginary parts of the eigenvalues A for the attractive 
case. The splitting of the lowest eigenvalues A — 4 and A = 16 
of the linear problem is highlighted by the vertical dashed line. 
The nonzero imaginary part in panel (c) corresponds to the 
eigenvalues originated by A = 16 of the linear problem for 
uj* < uj < 5 (see the text). Plots (d)-(f) are analogous to 
(a)-(c) but for the repulsive case. The branch studied is 
and the eigenvalues originated by A = 16 remain complex for 
the whole interval of ui studied (see panel (f)), and therefore 
the mode is unstable. 



C. Nonlinear modes of arbitrary amplitudes 

In order to study the linear stability of the nonlinear 
modes of a finite amplitude we have first calculated the 
modes using a modified shooting method developed in 
Ref. [llj and used in Refs. [HI, 0] to compute different 
families of nonlinear modes. Then we have found numer- 
ically the eigenvalues of the operator L~L+. To do so we 
have approximated the solution tp n (x) on a grid, replac- 
ing the second derivatives in by second-order finite 



7 




differences and caicuiated the eigenvalues of the result- 
ing sparse matrix. The result is shown in Figs. [5] and [31 
One can see that the families T^ 1 '^ at u = 5 possess two 
double eigenvalues (Fig. [21 panels b,c,e,f), A = 4 (merged 
eigenvalues A 2 ,i = A 2j3 ) and A = 16, (merged eigenval- 
ues A 2 ,o = A24). These eigenvalues split; in the case of 
A = 4 the resulting eigenvalues remain real and positive, 
one of them corresponding to the exact solution ([24)) . In 
the case of A = 16, if the nonlinearity is attractive, there 
exist two complex eigenvalues on the interval oj* < lu < 5, 
lu* Rj 3.83. At lu — lu* these eigenvalues merge and the 
mode becomes stable. Since the number of particles N of 
the mode grows when moving along the branch r 2 , one 
can say that there is a threshold on number of particles 
for the stability of the mode in attractive case. If the 
nonlinearity is repulsive, the two complex eigenvalues do 
not disappear through all the region of parameter lu in- 
vestigated; therefore the mode remains unstable. In the 
case of the family T 3 , (Fig. (a)-(f)) at lu — 7 there are 
three double eigenvalues A = 4, A — 16 and A = 36 which 
split. Then the scenario is similar to the case of the fam- 
ily r 2 . The eigenvalues originated by A = 4 remain real 
and one of them corresponds to the exact solution ([2"4"]l . 
In the attractive case the mode becomes stable after the 
two pairs of eigenvalues merge i.e. for u> < —2.10; at this 
point the second pair of the eigenvalues (both originated 
from A = 36) merges. In the repulsive case the mode 
remains unstable through all the region of the parameter 
lu studied. 

Summarizing, our results support that high- order 
modes of GPE with harmonic potential and attractive in- 



FIG. 4: Branches F„ and plots of the real and imaginary 
parts of eigenvalues A of the operator L^L„ vs to for V(x) = 
x 2 + 0.01a; 4 and n — 1. Plots (a)-(f) are organized as in 
Fig. [5] In the attractive case the dashed line on the plots 
(a)-(c) marks the upper boundary of the instability window. 
The lower boundary of this instability window is very close to 
N = 0, and not visible on the scale of the plots (a)-(f). Plots 
A, B and C show the branch T^ a ' and the real and imaginary 
parts of eigenvalues Ai,o and Ai,2 vs u) close to the linear limit 
uj = tii (i.e. N = 0) with magnification. On the plots A, B 
and C the lower boundary of the instability window is marked 
with a dashed line. 

teractions, are stable when the number of particles ex- 
ceeds a threshold value (different for each branch), what 
corroborates our analysis on the quasi-linear behavior of 
upper modes made in the begining of Sec. Ill Cl In the re- 
pulsive case, high-order modes of GPE with a harmonic 
potential are unstable. Our result contradicts those of 
|13| where it was claimed that stable modes exist for both 
signs of nonlinear term. 

IV. ANHARMONIC POTENTIALS (I) : SMALL 
PERTURBATIONS OF A HARMONIC 
POTENTIAL 

Now we turn our attention to the GPE with a harmonic 
potential perturbed by a quartic term, V(x) — x 2 + nx 4 , 
< I ft I -C 1. In this case the eigenvalues u) n and eigen- 
functions ip n {x), n = 0, 1, . . ., for the linear problem ([7]) 
can be found numerically or by means of asymptotic pro- 
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FIG. 5: Branches F„ and plots of the real and imaginary 
parts of eigenvalues A of the operator L~L^ vs lu for V(x) = 
x 2 + 0.01a; 4 and n = 2. In attractive case the lower boundary 
is very close to N = 0, and not visible on the scale of the 
figure. All the plots are organized as in Fig. [2] 




FIG. 6: Branches F„ and plots of the real and imaginary 
parts of eigenvalues A of the operator L~L^ vs lu for V(x) = 
x 2 + 0.01a; 4 and n = 3. In attractive case the lower boundary 
is very close to N = 0, and not visible on the scale of the 
figure. All the plots are organized as in Fig. [2] 



cedures [281 ]. Also one have to employ numerics for the 
construction of the branches of nonlinear modes and 

(r) 

r„ ; which are similar to the corresponding branches for 
harmonic potential case except for small deformations. 
As in the purely harmonic case, the branches r! a ' r ' ) are 
monotonic (at least for moderate values of N, lu and n) 
and can be parametrized by any of the parameters N and 
lu. However, the spectrum uu n is no longer equidistant, 
which has several important consequences as we will dis- 
cuss in what follows. 



A. Small amplitude modes 

The stability of small amplitude modes belonging to 
the branches T^ 1 '^ in the case of the anharmonic poten- 
tial can be also studied by means of asymptotic expan- 
sions (|T%l) - (|2"2"j) . Consider the case e = 0. The spectrum 
of the operator £„ now consists of simple eigenvalues 
A n ,fc = (Qk — u n ) 2 , k = 0, 1, . . ., and the Krein signature 
of the eigenvalue is if n .fc = sign(n — k). There- 

fore, if n is small the spectrum of £^ contains n pairs 
of close eigenvalues of opposite Krein signatures, A ni /- 
and A„ 2n-fcj fe = 0, 1, ..,n — 1, which are, nevertheless, 
different. The spectrum contains also a zero eigenvalue 
and an infinite sequence of increasing positive eigenval- 
ues. Therefore, generally speaking, the spectrum K n ^, 
k = 0, 1, . . . does not contain multiple eigenvalues. 

When passing from e = 0to0<e^;l the eigenvalues 
A„,fc of the operator L~L+ vary continuously. There- 
fore, one can expect the stability of the mode until two 
eigenvalues with opposite Krein signature merge. So, the 
small amplitude modes are expected to be stable for all 
the branches and rl r ' . 



B. Nonlinear modes of arbitrary amplitude: re > 0. 

In order to study the linear stability of the nonlinear 
modes of a finite amplitude we have first calculated the 
spectrum of the operator L~L+ numerically, concentrat- 
ing on the branches T^"'^ for n = 0, 1, 4. The obtained 
general picture appears to be very different from that of 
the harmonic potential. The plots of real and imaginary 
parts of the eigenvalues A vs lu for the branches Ti, T2 
and T3 in both attractive and repulsive cases are shown 
in Figs. H and© 

The general property of all the cases considered is that 
the instability of the nonlinear modes occurs only due to 
collisions of pairs of eigenvalues which are continuation of 
the respective eigenvalues A n> k and A„ : 2n-fe in the linear 
limit, i.e. of the eigenvalues of the operator C^. In what 
follows we call them (A„ fe, A nt 2 n —k)-pair8. The eigenval- 
ues in these pairs have opposite Krein signatures. 

Then, the other numerical results can be structured as 
follows. 
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a. Repulsive nonlinearity. The ground state modes 

(r) 

corresponding to the branch Tq are stable, since there 
are no (A^fc, A n ,2n-fc)-pair in the spectrum of LqLq. 

(r) 

The modes from the next branch, T\ , in the limit of 
strong nonlinearity correspond to "dark" soliton modes; 
there is one (A n< k, A„ ) 2 n -/s)~P au " m the spectrum of 
L^L^, but no collisions of eigenvalues have been found 
when tracing the modes of this branch within the range 
of parameter N where it has been investigated (see Fig. 
21 right panels). Therefore we conclude that the modes 
from are stable. A similar situation takes place for 

the branch r 2 where two (A„.fc, A„.2n-/c)-pairs in the 
spectrum of L 2 L 2 present: no collision has been ob- 
served within the range of parameter N under consider- 
ation (see Fig. [5j right panels). However this is not the 
case for higher branches. For instance, the collisions of 
eigenvalues has been observed for nonlinear modes from 

(r) — 4- 

I3 . The spectrum of the operator L 3 L3 includes three 
(A„ fc, A„ 2n-fc) _ P a i r and a collision of one of them (high- 
est) at some large enough value of N has been seen (see 
Fig- El right panels). After the point of collision (i.e. for 
greater values of ./V) the pair of collided eigenvalues be- 
come complex which means the instability of correspond- 
ing nonlinear modes. A similar situation takes place for 

(r) 

the branch T\ ' . In general, this points out to the fact 
that the instability of higher modes, generically, takes 
place, if the number of particles N exceed some thresh- 

(r) 

old value, which is particular for each branch r„ . The 

(r) 

existence of the threshold value for the branches T\ ' and 

(r) 

T 2 which we have not found in our numerical investiga- 
tion needs more delicate analysis. 

b. Attractive nonlinearity. 



The ground state 
modes of the branch Iq are stable. For the branch 

r^ Q) , there is one (A^, A n! 2n-fc)-pair in the spectrum of 
L\L\ (see Fig. [H left panels). Then there exist two 
bifurcation values of the parameter N (the number of 
particles). When N, increasing, reaches the first bifurca- 
tion value, the eigenvalues of this pair collide and become 
complex. For the values of N below this threshold, the 
mode is stable; however, for the potential under consid- 
eration the first bifurcation value is very tiny, so it is not 
visible in Fig. 2J Then, when N reaches the second bi- 
furcation value, the complex eigenvalue collide again and 
become real. So, we can conclude that the modes of 
are stable unless the number of particles N belongs to 
an instability window; the lower bound of this window is 
close to zero (but separated from zero). The size of the 
window of instability increases when k grows and both 
bounds of this window are quite sensitive to variation 
of k. A similar situation has been observed for higher 
branches of nonlinear modes. For instance, the spectrum 
of the operator spectrum of L 2 L 2l the branch I2 , in- 
cludes two (A n) fc, A nj2n -fc)-pairs. As N grows both of 
them undergo the same evolution: they become complex 
at first bifurcation value of N and return to be real at 
the second bifurcation value (see Fig. [5J left panels). 




FIG. 7: Branches rl a,r ^ and plots of the real and imaginary 
parts of eigenvalues A of the operator L~L^ vs lu for V(x) = 
x 2 — 0.01a; 4 and n = 1. All the plots are organized as in Fig. [2] 




FIG. 8: Branches ri?'^ and the plots of real and imaginary 
parts of eigenvalues A of the operator L^L„ vs lu for V(x) = 
x 4 and n — 1. All the plots are organized in the same manner 
as in the Fig. [2] 



The interval with respect to N between the first bifurca- 
tion value for the pair (Aa.i, ^-2,3) and second bifurcation 
value for the pair (A2 j o,A2 i 4) represents the window of 
instability. Upper boundary of this instability window is 
marked by dashed line in Fig. However, since the first 
bifurcation value for the pair ^2,1^2,3) (lowest curve 
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in Fig. [5l panel (b)) is very tiny, so the lower bound- 
ary of instability window cannot be separate from zero 
in Fig. [5] Therefore, the modes from are stable if 
N does not belong to instability window. This situation, 
probably, is generic for other higher branches. 

To confirm our results on the stability of nonlinear 
high-order modes we also have performed a series of di- 
rect numerical simulations of their evolution, perturbed 
by a random perturbation of 5% amplitude of the mode. 
Thus we have simulated the evolution of initial data of 
the form ipo(x) = ipn{x) [1 + r(x)] with r(x) a white noise 
of maximum amplitude 0.05. The subsequent dynamics 
of the modes under Eq. ([5]) was computed using a sec- 
ond order in time split-step pseudospectral scheme dis- 
cretized in space using trigonometric polynomials (via 
the FFT). In all the cases studied, which included most 
of the branches presented here our test verified the pre- 
dictions based on linear stability analysis. 



FIG. 9: Branches F„ and the plots of real and imaginary 
parts of eigenvalues A of the operator L~L^ vs lu for V(x) = 
x 4 and n — 2. All the plots are organized in the same manner 
as in the Fig. 





FIG. 10: Branches T„ and the plots of real and imaginary 
parts of eigenvalues A of the operator L^L^ vs lu for V(x) = 
x 4 and n — 3. All the plots are organized in the same manner 
as in the Fig. [2] The bold circle on the plot (e) highlights the 
collision of the eigenvalues of opposite Krein signatures which 
does not lead to instability. 



C. Nonlinear modes of arbitrary amplitude: ft < 0. 

Let us briefly summarize the stability results for the 
potential V(x) = x 2 + kx 4 , < |k| <C 1 and n < 0. 
In this case only a finite number of branches T^ 1 '^ can 
exist, since there is a finite number of discrete eigenvalues 
for Eq.([7]). These branches can be found numerically. 
The linear stability analysis performed for the potential 
V{x) = x 2 — O.Olir 4 shows (see Fig. [7]) that all solutions 

of the branch which we have considered are linearly 

(r) 

stable. On the other hand, the modes of T\ are also 
stable, except some instability window situated close to 
the point of branching. This is in contrast with the case 
k > 0, since in that case the instability window is situated 

(a) (r) 

on the branch T\ but not T\ . The solutions of the 
branch are also linearly stable whereas the solutions 

(r) 

of the branch T 2 are stable only in the vicinity of the 
branching point, i.e. only for JV < 1. For the branches 
T^' r ' the picture of stability/instability becomes more 
complex. 



V. ANHARMONIC POTENTIALS (II): 
POTENTIAL V[x) = x 4 

In order to show that the results for anharmonic po- 
tentials are quite generic, let us consider GPE with 
V(x) — x 4 . Again, the eigenvalues, Q n , and eigenfunc- 
tions tp n (x), n = 0,1,..., for the linear problem ([7]) 
cannot be obtained exactly. The branches of nonlinear 

(a) (r) 

modes T n and r„ which have been found numerically 
for n = 0, 1, 2, 3 are similar to the corresponding branches 
for harmonic potential case. These branches are mono- 
tonic (at least for moderate values of N, to and n) and 
can be parametrized by any of the parameters N , u>. 
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Due to the reasons discussed above, the small ampli- 
tude solutions of GPE with the potential V(x) = x A 
generically are stable for both, attractive and repulsive 
nonlinearities. In order to study the stability of nonlin- 
ear modes in general, we have analyzed the bifurcations of 
eigenvalues of L~L^ when ip n {x) varies along the families 
r„ and Tn\ Our numerical investigation shows that, 
as in the case of weak anharmonicity, the instability of 
nonlinear modes from T n occurs only due to collisions 
of those pairs of eigenvalues which are continuation of 
eigenvalues A n ^ and A„ 2n -fc in linear limit i.e. in the 
spectrum of C 2 n . These eigenvalues A n> k and A„^n-fc are 
of opposite Krein signature. Again, let us refer to these 
pairs of eigenvalues as (A n! fc, A n! 2n-/c)-pairs. 

Qualitatively, the stability picture for the case V(x) = 
x 4 is the same as in the case V(x) — x 2 + nx 4 , n > 0, 
except for some minor differences. The first one is that 
we have not found the threshold of instability in N in 
the case of repulsive nonlinearity for the modes in the 

(r) 

branches r„ , n = 0, 1, 2, 3. Therefore, these modes are 
stable in all the parameter range studied. Second, we 
have not found an upper bound of instability window in 



N in the case of attractive nonlinearity for the modes 
from r„ , n = 1,2,3. We believe that it is a technical 
problem related to the finiteness of the region studied 
but we think an upper bound of this instability window 
exists. The plots of real and imaginary parts of eigenval- 
ues A of the operator L~L+ versus ui for the branches 
Ti, T2 and for attractive and repulsive nonlinearities 
are shown in Figs. [H[9l andfTUl 

It is interesting to mention that we have also observed 
collisions of eigenvalues with opposite Krein signatures 
which belong to different (A n! fc, A rii 2n-fc)-pairs. In our 
case they did not lead to instability since they remained 
real after the collision (see for example Fig. [101 panel 
(e)). This phenomenon does not correspond to a generic 
situation; it is caused by the opposite parity of colliding 
eigenf unctions (one of them was odd, while the other one 
was even). 

We have also performed a set of direct numerical sim- 
ulations of the evolution of perturbed stationary modes 
in Eq. ([5]) as described in Sec. HVBl The outcome of 
those simulations confirms the results of the linear sta- 
bility analysis. 





V(x) = x 2 


V(x) = x 2 + KX 4 , ft > 


V{x) = x 4 


0-th mode, (ground 
state), a = ±1 


stable 


stable 


stable 


1-st mode, a — 1 small 
amplitude limit 


stable 


stable 


stable 


1-st mode ("bright soli- 
ton" ), a = 1 general case 


stable 


unstable, if iV belongs to some "in- 
stability window" . The lower bound of 
this "window" is separated from zero. 
The size of the "window" grows with k. 
Otherwise stable. 


unstable, if N belongs to some 
"instability window", with lower 
bound separated from zero. Oth- 
erwise stable. 


1-st mode, a = — 1 small 
amplitude limit 


stable 


stable 


stable 


1-st mode, ("dark soli- 
ton") a = —1, general 
case 


stable 


stable for all N which have been con- 
sidered. 


stable for all N which have been 
considered. 


Higher modes, a — 1 
small amplitude limit 


unstable 


stable 


stable 


Higher modes, a = 1, 
general case 


stable, if 

N exceeds a 
threshold. 


unstable, if JV belongs to some "insta- 
bility window", with lower bound sep- 
arated from zero. The size of the "win- 
dow" grows with k. Otherwise stable. 


Hypothetically, unstable, if N lies 
in a large "window" of instability. 
The upper boundary of this "win- 
dow" has not been found in our nu- 
merics. 


Higher modes, a = — 1 
small amplitude limit 


unstable 


stable 


stable 


Higher modes, a = — 1, 
general case 


unstable 


Hypothetically unstable, if N exceeds 
some threshold. Otherwise stable. 


stable for all N which we have con- 
sidered. 



TABLE III: Comparison of the stability properties for the potentials V(x) = x 2 , V{x) — x 2 + kx 4 (for k > 0) and V(x) — x 4 . 
The results for k < are not included in this table but discussed in the text 



VI. CONCLUSIONS AND DISCUSSION limit, the nonlinear WKB approximation, the Krein sig- 



Using a combination of different analytical and numer- 
ical tools including the analysis of the small amplitude 
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nature and direct numerical simulations we have ana- 
lyzed the stability properties of higher-order nonlinear 
trapped modes for the GPE with different potentials. 
First, we have reviewed the results for the harmonic po- 
tential V(x) — x 2 and discussed how the stability of the 
modes is essentially affected by the fact that levels are 
equidistant. Next, we have considered the weakly anhar- 
monic potential V(x) — x 2 + nx 4 , < \n\ <C 1. Our 
results, summarized in Table Hill lead to the conclusion 
that even a small anharmonicity which does not affect 
essentially the shape of the modes, improves drastically 
the stability properties of higher-order modes due to the 
fact that none of these potentials has an equidistant spec- 
trum. We conjecture that the same situation would take 
place also for more generic perturbation of harmonic po- 
tential, for instance, by non-symmetric (e.g. cubic) per- 
turbation. 

Then we have checked that in the case of stronger 
anharmonicity V(x) = x 4 the stability/instability pic- 
ture is similar to the case of potential V(x) — x 2 + nx 4 , 
< k <C 1, K > 0. We have studied the GPE with the 
potential V(x) = x 6 (the details have not been discussed 
in this paper) and found that they reproduce the same 
essential features. 

It follows form the arguments presented, that the sce- 
nario for appearance of instability induced by the equidis- 
tant spectrum of the harmonic oscillator holds also for 
other classes of potentials with equidistan t sp ectra (for 
construction of such potentials see [H [H, 34]) or, more 
generally, for potentials for which the spacing between 
some levels (not necessarily adjacent) are equal. In that 
situation, the splitting of double eigenvalues for the op- 
erator £^ can lead to complex eigenvalues in the linear 
stability problem. 

An interesting point for further investigation, is the ef- 
fect of the type of confining potential on the stability of 
higher order modes in two spatial dimensions, e.g. the 
stability of vortices under deformations of the potential. 
This subject has attracted a lot of attention in the last 



years [H GJ, El EI EI EI EH El S3 and the methodol- 
ogy developed in this paper could be useful. In fact, the 
situation is similar to the one considered above. In the 
case of harmonic potentials the spectrum of correspond- 
ing eigenvalue problem is equidistant; the correspond- 
ing eigenfunctions are Gauss-Laguerre modes. Then, one 
can expect that switching to anharmonic potentials can 
also change the stability properties of vortices and other 
higher order modes. 

Finally, we would like to mention another practical im- 
plication of the enhanced stability of nonlinear modes by 
the anharmonicity of the trap potential. As it was sug- 
gested in [n| such modes can grow from the eigenstates of 
the linear oscillator by increasing the nonlincarity using 
Feshbach resonance management (in the language of this 
paper this corresponds to the "motion" along a nonlinear 
branch starting from the bifurcation point as the number 
of particles increases starting from zero). This fact can be 
used for the generation of single solitons or even solitonic 
trains. The instability of the nonlinear modes in the case 
of the harmonic potential was the major obstacle for the 
practical implementation of that mechanism. However, 
the idea becomes experimentally feasible if an anharmoic 
potential is used since now higher order branches have a 
different stability and thus can lead to stable solitons. 
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